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Abstract 


The  AFIT  Total  Variation  Diminishing  Euler  Code  (ATEC)  was 
modified  to  include  a  thermal  nonequilibrium  model  to  investigate  high- 
temperature  effects  associated  with  vibrational  relaxation  in  a 
transonic  turbine  cascade.  Incorporation  of  this  model  into  ATEC  and 
creating  ANTEC  (AFIT  Nonequilibrium  TVD  Euler  Code)  was  accomplished  in 
three  phases.  Steady-state  solutions  obtained  with  ANTEC  were  compared 
with  those  obtained  with  ATEC  for  various  inlet  and  exit  conditions. 
The  CFL  criterion  was  held  constant  in  ATEC;  however,  it  required 
variation  for  ANTEC.  Blade  temperature  profiles,  temperature  difference 
contours  in  the  vicinity  of  the  trailing  edge  and  the  value  of  y  along 
the  blade  were  analyzed.  Even  when  corrected  for  high  temperatures,  the 
assumptions  of  a  calorically  perfect  gas  and  thus  a  constant  value  of  y 
are  inaccurate  due  to  the  temperature  dependent  nature  of,  Cp  and  Cy. 
Maximum  temperature  differences  of  -741°K  and  539°K  were  found  near  the 
trailing  edge  for  the  highest  temperature  case,  with  differences  being 
most  noticeable  through  the  expansion  at  the  trailing  edge  on  the 
pressure  surface  and  across  the  shocks  from  both  surfaces.  The 
vibrational  relaxation  model  showed  limitations  at  low  temperatures. 
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THERMAL  NONEQUILIBRIUM  EFFECTS  ON 


TURBINE  CASCADE  AERODYNAMICS 


_ Introduction 


I  ■  1  Background 

Current  algorithms  used  to  model  flow  through  transonic  turbine 
engine  cascades,  including  ATEC  (AFIT  Total  Variation  Diminishing  (TVD) 
Euler  Code),  assume  a  perfect  gas  and  neglect  high-temperature  effects 
such  as  molecular  vibration,  dissociation,  and  ionization.  Across  a 
shock  wave,  kinetic  eneray  is  converted  to  thermal  energy.  At  high 
temperatures  this  thermal  energy  can  vibrationally  excite  the  air 
molecules.  Ignoring  such  real  gas  effects  can  lead  to  an  extreme 
miscalculation  of  flow  temperatures.  For  example,  using  perfect  gas 
theory  for  a  blunt-nose  insulated  vehicle  at  60,  000ft  (T_,=  390°R/472°K) 

with  a  Mach  number  of  14,  a  temperature  of  over  10,000°K  is  predicted 
behind  the  shock  (the  surface  of  the  sun  is  6,  000°K)  .  The  temperature 
will  actually  be  4,900°K  because  molecules  absorb  energy  in  a  real  gas. 
(13:549)  Although  existing  transonic  turbine  cascades  do  not  run  at 
similar  flow  conditions,  flows  do  reach  temperatures  w.here  oxygen  begins 
to  dissociate  (2000-2500°K  at  one  atmosphere  of  pressure)  (2:19,374); 
above  500°K  the  specific  heat  at  constant  pressure,  Cp  ,  can  no  longer  be 
assumed  constant  for  diatomic  molecules  such  as  0^  and  N2 .  (3:63) 

Assuming  that  turbine  cascades  of  the  future  will  be  running  at  even 
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higher  temperatures  than  present  ones,  incorporation  of  a  thermal 
nonequilibrium  gas  model  into  ATEC  is  an  important  first  step  towards 
more  accurate  simulations  of  transonic  turbine  flow  fields. 

This  effort  is  a  merging  of  research  and  code  development  recently 
completed  at  AFIT.  ATEC,  developed  by  Capt  Mark  Driver,  is  'Vn 
explicit,  time-accurate,  two-dimensional,  finite-volume,  Euler  solver 
. . .  with  the  capability  of  resolving  the  complex  shock  structure  in  a 
typical  transonic  rotor  cascade."  (7:1)  The  code  uses  the  TVD  scheme 
developed  by  I'arten  as  modified  by  Yee.  (7,25)  ATEC  has  demonstrated 
improved  resolution  over  other  second-order  shock  capturing  schemes  such 
as  Lax-Wendrof f ,  and  comparisons  to  available  experimental  and 
analytical  solutions  have  been  favorable.  (7:2)  The  selected  thermal 
nonequilibrium  model  is  identical  to  that  developed  by  Capt  Ken  Moran  in 
h’  -  i  '  of  :  flov  ^bont  a  simple  axisymetric  blunt-body 
at  moderate  hypersonic  "peeds.  This  mode]  was  developed  for  possible 
use  with  mu’ t i-species  gases.  (16,17)  The  purpose  of  thxs  study  is  to 
investigate  the  thermal  nonequilibrium  effects  and  compare  results 
assuming  a  perfect  gas  for  turbine  carcaoes. 

1^ _ Outline  of  Study 

As  with  any  problem,  breaking  this  research  into  smaller  pieces 
facilitated  progress  and  understanding.  Phase  One  consisted  of 
implementing  a  second-order  TVD  scheme  for  a  shock  tube  problem. 
Harten's  ULTIC  algorithm  (12)  was  used  due  to  its  capability  to 
accurately  capture  the  shock  and  the  contact  surface.  Harten  also 
presents  shock  tube  data  with  which  comparisons  could  be  made. 
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Although  there  is  a  strong  interaction  between  molecular  energy 
modes  (15:4),  a  simpler  model  was  investigated  first;  therefore,  only 
the  vibrational  relaxation  model  was  incorporated  into  the  shoc)c  tube 
problem  for  Phase  Two.  Moran  investigated  two  vibrational  relaxation 
models:  the  Milliken  and  White  model  and  the  Park  model.  (16:22-24) 
The  Milliken  and  White  model  computes  relaxation  times  based  on  an 
empirical  formula  (Equations  9  and  10) .  The  Park  model  incorporates  the 
collision  time  and  was  developed  for  suborbital  hypersonic  flight. 
(19:440)  The  applicability  of  the  Park  model  at  moderate  to  low  Mach 
numbers  is  not  known.  Therefore,  the  Milliken  and  White  model  was 
chosen  for  this  investigation. 

The  TVD  scheme  used  for  the  first  and  seco-.^  phases,  ULTIC,  is 
based  on  a  second-order  accurate,  upwind-dif ferencina  scheme.  The 
unique  character  of  TVD  schemes  is  the  numerical  dissipation  in  the  flux 
term.  According  to  Yee  (25),  TVD  schemes  were  developed  for  homogeneous 
hyperbolic  problems  and  a  nonequilibrium  source  term  may  or  may  not 
preserve  the  original  TVD  properties.  However,  if  each  species  of  a  gas 
is  assumed  to  behave  as  a  thermally  perfect  gas,  the  flux  function  will 
possess  the  homogeneous  property.  (25:117)  It  is  this  flux  function 
that  makes  the  scheme  stable  and  robust  near  flow  field  discontinuities, 
resulting  in  good  shock  capturing  capability  The  dissipation  is 
locally  adjusted  and  applied  differently  to  the  linearly  degenerate 
fields  and  the  genuinely  nonlinear  ones.  Linearly  degenerate  fields  are 
those  in  which  =  0  and  are  exclusively  contact  discontinuities; 
where  a*'  represents  an  eigenvalue,  R*'  represents  an  eigenvector  of  the 
Jacobian  matrix  (explained  in  Section  II. 2),  and  k  is  valued  from  1  to  5 
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in  this  study. 


k  k 

In  contrast,  nonlinear  fields  (a  R  ^0)  are  either 
shocks  or  rarefaction  waves.  (12:374)  With  this  form  of  numerical 
dissipation,  the  smearing  of  the  contact  surface  and  leading  and 
trailing  edges  of  the  expansion  fan  is  not  as  pronounced  as  with  a 
first-order  Roe  scheme  or  the  other  second-order  TVD  schemes  discussed 
by  Harten.  (i2: 380-393)  The  linearly  degenerate  fields  are  modified  to 
be  slightly  convergent,  thus  reducing  smearing  and  loss  of  resolution  in 
the  contact  discontinuity  region.  The  characteristic  fields 
corresponding  to  the  eigenvalue  u  are  linearly  degenerate  and  the  fields 
corresponding  to  u ± C  are  genuinely  nonlinear.  (12:379) 

Phase  Three  saw  the  merging  of  the  relaxation  model  into  ATEC, 
which  uses  Yee ' s  modification  to  ULTIC  due  to  the  over  diffusive  nature 
of  ULTIC  in  two  dimensions.  (8,25:25)  This  is  also  a  shock  capturing 
scheme  that  eliminates  the  need  to  know  apriori  the  complex  shock 
structure  in  a  turbine  cascade.  In  general,  classical  shock  capturing 
schemes  have  resulted  in  non-physical  oscillations  near  discontinuities. 
In  contrast,  these  oscillations  are  not  observed  with  TVD  schemes.  The 
modified  version  will  be  referred  to  as  ANTEC  (AFIT  Nonequilibrium  TVD 
Euler  Code) . 

At  this  point,  comparisons  between  ATEC  and  ANTEC  were  made  at 
various  flow  conditions  to  determine  any  'high-temperature'  effects. 
First,  a  low-temperature,  low-pressure  case  was  used  for  debugging  and 
verification  of  ANTEC.  Even  though  the  blade  geometry  was  for  a 
transonic  rotor  cascade  designed  by  NASA  and  not  that  of  a  FlOO  turbine 
blade,  the  FIOO-PW  turbine's  inlet  temperature,  pressure,  and  stage 
pressure  ratio  were  chosen  for  the  next  case  as  being  representative  of 
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turbine  conditions. 


Actual  numbers  for  turbines  currently  in  design 


were  not  available  due  to  their  sensitive  nature.  This  research  was 
concerned  with  the  concept  of  thermal  nonequilibrium  effects  in  the  flow 
through  a  turbine  cascade.  The  NASA  blade  geometry  and  the  FIOO-PW 
turbine  conditions  simulate  a  hypothetical  turbine  for  research 
purposes.  Once  the  FlOO  inlet  case  was  run,  the  inlet  temperature  was 
increased,  holding  the  other  factors  constant  to  investigate  any  high- 
temperature  effects;  such  as  temperature  differences  between  the  two 
codes.  These  differences  were  seen  to  be  greatest  in  the  vicinity  of 
the  trailing  edge. 

_ Assumptions 

Throughout  this  study  certain  assumptions  were  made  to  simplify 
the  problem; 

1.  Air  is  a  single  "composite"  species  approximately  consisting  of 
21%  O2  and  79%  Nj  with  a  molecular  weight  of  28.9844  kg/kg- 
mole.  (20:87) 

2.  The  gas  is  thermally  perfect. 

3.  The  effects  of  molecular  transport  (viscosity,  diffusion, 
thermal  conductivity)  are  negligible.  (22:178) 

4.  There  are  no  dissociation,  ionization  or  electronic  effects. 

5.  All  energy  modes  decouple. 

For  an  analysis  based  on  a  multi-species  gas,  a  continuity 
equation  and  vibrational  energy  equation  are  necessary  for  each  species. 
(16,17)  Therefore,  the  first  assumption  simplifies  the  structure  of  the 
governing  equations.  A  thermally  perfect  gas  has  a  constant  value  of  R; 
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however,  the  specific  heats,  Cp  and  Cy,  and  their  ratio,  7,  are  not 

necessarily  constant.  The  assumption  of  a  thermally  perfect  gas  also 
allows  retention  of  the  desired  TVD  properties.  (25:117)  The  remaining 
assumptions  serve  to  more  fully  isolate  nonequilibrium  effects  due  to 
molecular  vibration.  The  degree  of  dissociation  of  a  gas  will  decrease 
as  pressure  increases,  making  the  fourth  assumption  valid  for  a  high- 
temperature,  high-pressure  turbine  cascade  (approximately  20  atmospheres 
at  the  inlet).  (22:152-157) 

Yee  discusses  second-order  TVD  applications  for  nonequilibriurr; 
flows  (25:115-117)  which  were  applied  by  Moran  in  his  research  of 
thermal  and  chemical  nonequilibrium  flows  in  a  one-dimensional  shock 
tube.  (17)  Unfortunately,  no  literature  was  found  that  dealt  with 
nonequilibrium  flow  through  a  turbine  cascade.  Also  limited 
experimental  data  for  flow  through  a  turbine  cascade  was  found  that 
could  possibly  be  used  for  comparisons.  (6,9)  In  the  data  that  was 
found,  there  were  significant  radial  pressure  gradients  across  the  duct. 
ATEC  could  be  modified  to  account  for  variations  at  the  inlet;  however, 
it  was  felt  that  a  mure  systematic  comparison  of  ATEC  and  ANTEC  could 
be  obtained  without  changing  the  inflow  conditions. 
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Since  a  weak  solution  is  anticipated,  the  conservative  or 


divergence  form  of  the  governing  Euler  equations  is  used.  This  choice 
also  simplifies  the  derivation  of  a  second-order  explicit  algorithm. 
(23:87)  The  governing  equations  in  conservation  form  are 


^  9  E(U)  9F(U) 

8 1  9  X  9  y 

where 


p 

"  pu 

pv 

"  0  ■ 

P  u 

p  u^+P 

p  uv 

0 

u  = 

pv 

E, 

E  = 

puv 

(E,+P)u 

F  = 

pv^+P 

(E,+P)v 

S  = 

0 

0 

w 

—  Evib  — 

—  Eyitjd  — 

—  Eyib'^  — 

These  are  the  Euler  equations  for  conservation  of  mass,  momentum 
and  energy  for  an  inviscid,  non-conducting  gas  plus  a  fifth  equation  to 
account  for  the  decoupled  vibrational  energy.  The  source  term, 

(defined  in  Equation  8) ,  is  the  Landau-Teller  form  describing  the  energy 
exchange  between  energy  modes,  which  assumes  that  molecules  are  free  to 
vibrate  in  any  quantum  level,  but  changes  in  energy  only  occur  with 


adjacent  levels.  (3:65,17:8)  E,  is  the  total  energy  per  unit  volume  as 
defined  by 


E, 


_p_,£ 

(ri)  2 


(U^+V^) 


Eyitj 


(3) 


is  the  vibrational  energy  per  unit  volume  (E^,^  =  .  The 
equilibrium  model  characterizing  vibrational  energy,  is  based  on 
statistical  and  quantum  mechanics.  (22:86-139)  This  model  states  that 
all  internal  energy  is  partitioned  into  translational,  rotational, 
vibrational  and  electronic  modes.  The  equilibrium  vibrational  energy  is 
given  by  the  equation 


eyib  = 


a  (In  Q^,t,) 


ai 


(4) 


with  ^vib  being  the  partition  function 


^vib 


1 


1 


-ev 

exp  — 


(5) 


After  substitution  of  Equation  5  into  the  equilibrium  equation  for 
vibrational  energy.  Equation  4  becomes 


R  0„ 


®vib  ~ 


exp 


(6) 
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The  characteristic  vibrational  temperature,  0^,,  is  a  constant 
obtained  from  spectroscopic  data.  The  value  used  in  this  study  (0^  = 

3154. 8°K)  is  a  mole  weighted  value  from  the  assumption  made  in  Section 
1.3  that  air  is  made  up  of  0^  (0^  =  2270°K)  and  (0^  =  3390°K)  .  (22:135) 

Wvib  is  a  function  of  the  relaxation  time,  T.  Relaxation  ti.me  can 
be  considered  the  time  necessary  for  molecules  in  nonequilibrium  (for 
example,  air  through  a  shock)  to  absorb  energy  and  reach  an  equilibrium 
state.  The  rate  at  which  the  gas  will  approach  equilibrium  is  directly 
proportional  to  its  departure  from  equilibrium.  (11:196) 

As  with  Moran's  previous  nonequilibrium  analyses  (16,17),  this 
Study  will  use  a  two-temperature  model.  The  first  temperature,  T, 
identifies  the  translational/rotational  energy  in  equilibrium  and  the 
second  temperature,  identifies  the  nonequilibrated  vibrational 
energy.  Equilibrium  exists  when  T  and  T^ij,  are  equal,  in  which  case  the 
source  term,  will  be  zero. 


'  vib 


0v 


R  0, 

ln(-—  +  1) 

6uih 


(7) 


W  is  defined  as  (17:8) 


=  P 


"'vibj 


(8) 


with  being  the  vibrational  energy  defined  by  Equation  6  and  is 

'  vib 

extracted  from  the  solution  of  Vibrational  relaxation  timies  have 
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been  determined  experimentally;  and  below  5000°K,  it  has  been  shown  that 


for  the  most  common  species  in  air,  i  can  be  determined  by  (18;  58) 


exp  (A  (T-'^-.015p'^^)  -  18.42) 
P 


(9) 


A  =  1.16x10'^ 


(10) 


is  called  the  reduced  mass  of  the  air.  For  the  single  species  case, 
it  becomes  one-half  of  the  molecular  weight.  T  defines  the 
translational/rotational  energy  and  is  given  by 


T 


(11) 


P  =  (y-1)[E,  -  -  I  (u^  +  v2)]  (12) 

In  determining  t,  P  is  pressure  in  atmospheres.  Pressures  in  the  codes 
were  in  N/m^,  so  that  a  conversion  factor  of  101,325  N/m^/atm  was 
implemented  to  incorporate  Equation  9  into  ANTEC .  Clearly,  as  P 
increases  the  relaxation  time  will  decrease. 

R  is  the  gas  constant  defined  by 


R  = 


_ 2 _ 

molecular  weight 


(13) 


X  being  the  Universal  Gas  Constant  equal  to  8134  J/°K  mole. 
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Remembering  that  the  value  of  Cp  is  temperature  dependant  and  with 

the  assumption  of  a  thermally  perfect  gas,  an  equilibrium  or  effective 
value  of  Y  can  be  calculated  using  the  following  relationships 

Y  =  (14) 

and 

Cp  =  c,  +  R  (15) 

where 


c„  =  c. 


+  c. 


(16) 


Y  can  be  solved  for  using 

=  ^  R 

^^Trans/Rol  2 

(22:124,133)  and 


'^Tvib 


=  R 


\  2  exp  (Gy/T) 

.TJ  (exp  (e^/T)  -  1)2 


(17b) 


(22:135)  In  the  limit  as  T  approaches  zero,  also  goes  to  zero 
resulting  in  y  approaching  1.4.  As  T  increases,  the  value  of  y 
decreases;  and  as  T  becomes  larger  than  approaches  R, 
indicating  a  fully  excited  state.  If  the  temperature  were  to  continue 
to  increase,  the  energy  would  have  to  go  into  alternate  modes  such  as 
electronic  or  dissociation.  In  equilibrium,  Cy  is  dependant  upon  only 
one  temperature,  T.  However,  for  this  research,  a  generalization  for 
application  to  nonequilibrium  flow  is  necessary  for  comparison  and  Cy^^ 
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(Equation  17)  will  be  defined  using  Ty,[,  (Equation  7)  .  With  these,  a 
value  of  Y  can  be  constructed  for  a  nonequilibrium  flow.  This  will  be 
referred  to  as  Yg^.  The  values  of  Y  actually  used  in  the  codes  will  be 
referred  to  as  Ya  ATEC  and  Yno  ANTEC.  For  equations  that  are 
applicable  to  both  codes,  Y  will  be  unsubscripted .  When  the  equation  is 
applied  for  ATEC,  Ya  will  be  used;  and  conversely,  when  applied  to 
ANTEC,  Yne  will  be  used. 

Again,  this  research  is  only  concerned  with  vibrational  energy  and 
relaxation.  However,  it  should  be  noted,  as  a  point  of  reference,  that 
translational  equilibrium  is  reached  within  a  few  collisions  and 
rotational  equilibrium  is  attained  in  less  than  ten  collisions  for  most 
polyatomic  gases;  but  vibrational  equilibrium  may  require  thousands  of 
collisions.  (14:25,3:62)  So,  for  practical  purposes,  translational  and 
rotational  states  can  be  considered  to  be  in  equilibrium 
instantaneously,  whereas,  vibrational  relaxation  requires  a  finite 
amount  of  time,  which  is  approximated  by  t. 

Appendix  B  briefly  summarizes  the  modifications  made  to  ATEC  to 
create  ANTEC. 

II. 2  Finite-Difference  Scheme 

TVD  is  a  class  of  modern  higher-order  shoc)c  capturing  algorithms 
which  uses  nonlinear  numerical  dissipation  to  yield  a  stable, 
nonoscillatory  solution.  (25)  Yee  emphatically  points  out  that  the  TVD 
property  is  only  valid  for  homogeneous  scalar  hyperbolic  problems. 
However,  the  TVD  method  has  been  successfully  extended  to  the  treatment 
of  nonlinear  systems  of  hyperbolic  conservation  laws.  (25,7,17) 
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In  the  formulation  of  a  TVD  scheme,  it  is  necessary  to  locally 
obtain  the  set  of  eigenvalues,  a,  the  eigenvector  matrix,  R,  and  its 
inverse,  R“',  for  the  two  Jacobian  matrices: 


^  ^E(U) 

au 


B 


3f(U) 

au 


{ ISa, b) 


Fortunately,  these  quantities  are  available  in  the  literature. 
(7,16,25,12,15)  However,  it  is  also  possible  to  solve  for  a,  R,  and  R“' 
using  Mathematics  if  necessary.  The  eigenstructure  of  the  original 
version  of  ATEC  was  extracted  from  Yee.  (25,7,8)  The  nonequilibrium, 
version  incorporated  the  eigenstructure  presented  by  Moran.  (16:79-80) 
The  eigenvalues  for  A,  following  a  general  coordinate  transformation, 
are 


a  = 


Uc- 

LUc  +  K^C, 


(19) 


The  eigenvector  matrix  of  A  is 


R 


0 

-k. 


0  1 

0  u, 

0  V. 


Uc-k,c 


Vc-kzC 


Uc+k,c 


Vc+kzC 


(k.Vc-kjU,)  1  H-k,ux-k,v-C  H+k,u,c+k,v,c 


0 


‘'vib 


(20) 
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and  its  inverse 


k^u+kjV 

-kr 

ki 

0 

0 

■®vib^l 

u(e,ibb,i 

'  -evibbj 

"•  +®vibl 

1-bi 

ubj 

vbj 

-b3 

b3 

r’  = 

k.u  k,v 

b,+— +— 

1  c  c 

-b3U-- 

-b3v-- 

b3 

-b3 

2 

.  2  . 

L  2  . 

2 

2 

k.u  k,v' 

^-T-T 

k,l 

-‘'3+7 

fh  y 

■^3+c 

b3 

-b3 

2 

2 

2 

2 

2 

(21) 


where 


^3 


(V  -  1 ) 


b,  = 


22,23) 


k.  = 


2.i 


^2  ~ 


-k. 


e;  +  4;)’ 


(24,25) 


The  contravariant  velocity  in  the  direction,  u^.,  is 

determined  by 


(26) 


(27) 


where  the  speed  of  sound  is  C,  defined  by 


Uc  = 


= 


O' 
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(7-1)  H  -  |(u2  +  v2)  - 


(28) 


= 

Enthalpy,  H,  is  given  by 


(E,  +  P) 
P 


(29) 


In  the  'n-coordinate  direction,  v^.  is  used  in  place  of  u^,  and  ri  replaces 
^  in  all  derivatives  (Equations  19  and  24-27)  .  The  eigenvector  matri;-: 
of  B  and  its  inverse  have  the  same  form  as  Equations  20  and  21,  again 
using  T)-coordinate  derivatives  in  Equations  24  and  25. 

With  the  solution  being  generated  at  cell  centers.  Roe  averaging 
was  applied  throughout  the  analysis  since  the  eigenvalues  and 
eigenvectors  require  some  form  of  averaging  at  cell  interfaces.  This  is 
the  most  common  form  of  averaging  due  to  its  simplicity  and  its  ability 
to  resolve  discontinuities.  (25:59,6:5)  Figure  1  depicts  the  cell 
interfaces  (i,j)  ,  cell  centers  (i+i/2,j+i/2)  ,  and  the  cell  centers  outside  of 
the  boundaries  )cnown  as  ghost  points.  Roe  averaging  in  one  dimension  is 
implemented  by 


+  y-i 

Vi'2  D  +  1 


(30) 


where  x  is  u,  v,  H  or  e  ,  and  D  is  defined  as  follows 


D 


I  P. 


(31) 
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Although  the  form  of  D  chosen  was  for  a  perfect  gas  in  equilibrium,  i 


can  still  be  used  in  this  application  due  to  the  assumption  of  a 
thermalry  perfect  gas.  (25:59)  This  method  of  averaging  was  also 
discussed  for  nonequilibrium  flows  by  Liu  and  Vinokur.  (15:17) 


0  ^ 

0  ^ 

0 

i.i+1 

■-*  i  't-I  ' 

0 

Cl 

o 

i+1/2 .1  +  1/2 

1. 

j 

l.i  ^ 

0 

O 

0 

Cell  Ceniers  CXitside  o(  Boundary  (Ghost  Points) 


Figure  1 :  Cell  Interfaces  and  Cell  Centers 


with  X  defined  as  A^/A^  and  A(-S  is  the  treatment  of  the  source  term. 

The  n  superscript  is  now  dropped  for  simplicity.  The  numerical  flux, 

f  ,  is  evaluated  at  cell  interfaces  and  is  expressed  as 
1*1/2 


f 

]*I/2 


f(u;.f(u,,)  r;:., 

K=1 


(33) 


For  this  research  m  was  5.  2  are  the  eigenvectors;  and  the 

effective  numerical  viscosity,  P  ,  is  defined  as 

jtl/2 


34) 


with 


^'2 


«  k  K 

^a,*i/2  +  Y,* 


1*1/2 


-v 

V,  =  X  3; 


‘j*I/2 


(35) 

(36) 


and 


®j*l/2  ~  [^1*1/2]  ^(*1/2'-'' 


(37) 


=  U,.,  -  U, 


(38) 


Q,  the  entropy  correction  factor  or  also  )cnown  as  the  coefficient  of 
numerical  viscosity  (12:363),  and  a  are  functions  expressed  as 


o''^,,;(Z)  =  ^(Qiz)  -z^) 


(39) 


Q(Z)  =  —  +  £  for  |z|  <  2e 


4e 


(40) 


Q(Z)  =  lz|  for  iz|  >  '£ 
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e is  a  variable  parameter  in  the  code  that  can  be  selected  by  the  user. 


Values  between  0  and  1/2  are  allowed;  with  0  yielding  the  least 
dissipative  results.  (12:375,25:16) 

k 

For  the  linearly  degenerate  fields  ,  the  limiter,  is  determined 

from 


g,'  =  sf 


j+l/2 


max[0,min{cii’‘^,,2ja^ 


Y\a\ 


■S).,/2aj.„2,min(|aj 

♦id 


(41) 


with 

Sj1i/2  =  signd,  ajl,,^)  (42) 

The  'sign'  function  assigns  the  sign  (+  or  -)  of  the  second  term  in  the 
parentheses  to  the  value  of  1. 

(0  is  also  an  input  parameter.  Input  into  the  codes  is  explained 
in  Appendix  C. 

For  the  nonlinear  fields 


“j  i/2  ^“♦1/:  “i-vzl 


-k 

g, 


®j*i/2  ®i-i;2 


(43) 


With  Equations  37,  39,  and  43,  obtained  from 


(^1  -  g,') 


.„2~ 


for  *  0 


■j4l/2 


{44a) 


tvi  -  0 


for  =  0 


(44b) 
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Here,  Y  is  not  related  to  the  ratio  of  specific  heats;  use  of  this 
symbology  is  made  to  remain  consistent  with  the  literature.  (25,12,7) 

The  Strang-type,  fractional  step  method  is  incorporated  to  advance 
the  solution  in  time  (7:5) 


,  ,n+2  0^  0*'  0^  o*'  0^  ,  ,n 

'-*('.))  " 

where  £  is  an  operator  representing  Equation  32,  h  =  and  ^  and  T) 
correspond  to  the  direction  of  the  computational  sweep.  With  this 


method,  each  iteration  encompasses  two  time  steps. 


The  computational  grid  for  this  study  was  a  C-type  with  170x20 
points.  Since  averaging  at  cell  centers  was  required,  a  grid  of  cell 
centers  including  'ghost  points'  below  the  blade  and  branch  cut  at  the 
trailing  edge  was  generated  by  the  codes.  All  data  output  in  this 
research  were  values  at  cell  centers.  (See  Figure  1)  The  blade  geometry 
used  for  this  research  was  the  rounded  trailing  edge  blade  used  by 
Driver  and  Beran  as  shown  in  Figure  2.  (7) 
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FIGURE  2:  Grid  for  Rounded  Trailing  Edge  Blade 

II. 4 _ Boundary  Conditions  (7:3-5) 

Both  codes  used  the  same  boundary  conditions;  although  an 

additional  boundary  condition  was  incorporated  to  specify  in  ANTEC . 

The  subsonic  inlet  was  assumed  to  be  part  of  an  imaginary  duct 
extending  infinitely  far  upstream.  Disturbances  in  the  computational 
domain  are  allowed  to  pass  through  the  inlet,  without  reflection,  and 
travel  upstream  as  a  simple  wave.  This  assumption  motivated  the  use  of 
one-dimensional  characteristic  wave  theory  to  establish  boundary 
conditions  at  the  inlet.  One  characteristic  runs  from  the  domain 
interior;  therefore,  only  one  quantity  may  be  extrapolated  while  all 
others  must  be  specified.  The  speed  of  sound  was  extrapolated  while  the 
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(inlet  Mach),  Pc,  and  were  specified.  Once  the  speed  of  sound  at 
the  inlet,  C2,  was  determined,  Tj  was  derived  from  the  Riemann  invariant 


(  7  -  1) 


(7  -  1 ) 


(46) 


Far  upstream  the  flow  was  at  rest,  and  Uc  vanishes.  After  a  minor 
amount  of  manipulation.  Equation  46  becomes 


h 

T. 


(47) 


Given 

M2 

and  Tc, 

the 

inlet 

temperature,  Tj,  can  be  solved 

for . 

Once 

this 

was 

known, 

the 

inlet 

pressure  was  easily  obtained 

using 

the 

isentropic  relationship 


Y-i 


(48) 


The  inlet  vibrational  energy  necessary  for  ANTEC  was  also  determined  by 
Tj  with  the  assumption  that  at  steady  state  all  waves  will  have  passed 
back  through  the  inlet,  leaving  the  inlet  in  equilibrium. 

Simple  wave  theory  was  again  used  to  solve  for  the  e.xit 
conditions.  The  flow  was  assumed  to  exit  into  a  plenum  that  required 
the  exit  pressure,  P3,  to  match  the  plenum  pressure;  any  pressure 
aisturbances  were  reflected  back  into  the  computational  domain.  With 
the  exit  pressure  specified,  density  at  the  exit  was  calculated  using 
isentropic  relationships 
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(49) 


P3  = 


■’mu 


and 


fP, 


®int 


_ int 

Pint. 


(50) 


Entropy,  tangential  velocity,  and  the  Riemann  invariant  were 
extrapolated  from  the  interior  of  the  domain  using  one  point  and 
vibrational  energy  was  extrapolated  using  two  points. 

For  the  C-grid  used  as  shown  in  Figure  2,  periodicity  was  applied 
along  the  outer  boundary  since  the  computational  domain  models  one  blade 
in  an  infinite  cascade  of  blades.  These  periodic  boundary  conditions 
were  applied  at  the  ghost  points  located  outside  the  computational 
domain.  Along  the  wake  cut,  continuity  was  enforced. 

On  the  blade  surface,  an  adiabatic  wall  was  assumed  and  normal- 
momentum  was  used  to  calculate  pressure  and  velocities.  To  solve  for 
the  vibrational  energy  at  the  blade  surface  for  ANTEC,  was 
reflected.  (4) 


II. 5 _ Initial  Conditions  (7:5) 

The  solutions  for  both  ATEC  and  ANTEC  were  started  at  rest  with 
the  initial  conditions  for  p,  P,  and  T  in  the  computational  domain 
being  set  to  the  constant  values  far  upstream  of  the  inlet.  Due  to  the 
gas  in  the  cavity  starting  at  zero  velocity  with  no  discontinuities, 
equilibrium  was  assumed  to  exist  initially  to  specify  E^,j,  in  ANTEC. 

For  this  research,  three  cases  were  investigated  and  are 
summarized  in  Table  1: 
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1.  a  low-temperature,  low-pressure  case  {as  in  reference  1) 


2.  the  FIOO-PW  inlet  conditions  and  pressure  ratio 

3.  increased  inlet  temperature  with  FIOO-PW  pressure 
specifications . 

In  Cases  2  and  3,  was  adjusted  based  on  the  desired  inlet 
temperature,  Tj,  and  =  1.4  for  Case  1;  even  though  at  this  low 

temperature,  a  corrected  value  of  1.3968  could  have  been  used.  (26:703- 
704)  Decoupling  of  the  vibrational  energy  results  in  a  constant  value 
of  1.4  for  ANTEC .  This  is  based  on  the  analysis  of  Liu  and  Vinokur. 
(15)  They  define  a  pressure  derivative,  k,  which  is  necessary  in 
simplifying  the  nonequilibrium  eigenstructure  presented  in  the 
literature.  (15,16)  As  electronic  effects  are  excluded  from  this  study, 
K  is  defined  only  in  terms  of  the  translational/rotational  value  of 
(Equation  17a) ;  and  nonequilibrium  internal  energies  do  not  contribute 
to  K.  (15:9)  Therefore,  K  is  equal  to  2/5  or  Yne  "  ^  when  Yne  =  ^  • 


TABLE  1 

CASE  SUMMARY  FOR  DESIRED  INLET  CONDITIONS 


Case 

P2  (N/m2) 

Tj  (°K) 

1 

19.81x10*’ 

375.73 

2 

215.81x10*’ 

1682.78 

3 

215.81x10*’ 

2500.00 

An  inlet  Mach  number,  Mj,  was  calculated  with  the  following 
relationship : 
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2 


(Y  +  1 ) 


(51) 


■_W' 

w 

l  a) 


2  +  (Y  -  1) 


where  the  mean  design  W/W^,,  for  the  blade  is  0.317.  (21)  is  the 
ratio  of  the  relative  gas  velocity  to  the  relative  gas  velocity 
corresponding  to  a  Mach  number  of  unity.  Table  2  shows  the  values  of  y 
for  each  case  and  the  resultant  Mj. 


TABLE  2 

CASE  SUMMARY  FOR  INLET  MACH  AND  y 


Case 

M2 

y 

1 

ATEC 

0.352 

1.4000 

ANTEC 

0.352 

1.4000 

2 

ATEC 

0.297 

1.3053 

ANTEC 

0.292 

1.4000 

3 

ATEC 

0.298 

1.2896 

ANTEC 

0.292 

1 .  4000 

Using  Equations  47  and  48,  the  values  of  and  P_  for  Cases  2  and 
3  were  obtained  using  desired  inlet  conditions.  The  appropriate  values 
of  Mj  and  y  are  shown  in  Table  2.  The  values  for  and  used  in 
Case  1  were  calculated  by  Driver  and  Beran  and  given  to  the  author  as  a 
test  case.  Using  W/Wj.^  =  .381  as  in  reference  7,  values  of  Mj,  T,  and 
P2  for  Case  1  were  calculated  using  Equations  47,  48,  51  and  the  gas 
tables  in  reference  26.  The  value  of  Tj  used  in  Case  3  was  2500°K.  The 
pressure  ratio,  P2/P3/  was  calculated  to  be  3.11  for  Case  1  and  the  FlOO 
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pressure  ratio  used  was  3.66, 


Since  the  values  of  7^  and  Mj  for  ANTEC 


do  not  change  for  Cases  2  and  3,  the  ratio,  was  constant . 
Therefore,  P../P2  was  also  constant.  Emphasis  was  placed  on  matching  the 
inlet  conditions  between  the  two  version  so  that  accurate  comparisons 
could  be  made  and  appropriate  conclusions  drawn.  The  conditions 
calculated  far  upstream  are  summarized  in  Table  3. 


TABLE  3 

CASE  SUMMARY  FOR  CONDITIONS  FAR  UPSTREAM 


Case 

P_  (N/m2) 

T„  ("K) 

1 

ATEC 

29.31x10'’ 

420.16 

ANTEC 

29.31x10'’ 

420.16 

2 

ATEC 

312.57x10'’ 

1835 . 07 

ANTEC 

320.24x10'’ 

1883.64 

3 

ATEC 

311.30x10'’ 

2714.37 

ANTEC 

320.24x10'’ 

2798.41 
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III.  Results 


A  brief  discussion  of  the  shock  tube  results  can  be  found  in 
Appendix  A.  This  section  will  focus  on  the  comparative  results  obtained 
with  ANTEC  and  ATEC . 

III.l  Case  Summary 

The  cases  discussco  in  Section  II. 5  were  investigated  at  steady- 
state,  since  reason'  j±e  comparisons  could  only  be  made  between  steady- 
state  solutioto.  Once  attained,  the  flow  parameters  could  then  be 
evaluated.  Of  particular  interest  were  the  flow  temperatures  along  the 
blade,  in  the  vicinity  of  the  shock  and  the  trailing  edge. 

Table  4  shows  the  steady-state  values  attained  for  each  case  and 
the  percentage  difference  from  the  desired  inlet  values.  Both  versions 
did  poorly  compared  to  the  calculated  in  Table  2.  The  results  were 
closer  to  the  FlOO  turbine  inlet  Mach  number  of  0.4.  is  driven  by 
the  pressure  ratio;  or  if  the  flow  is  choked,  it  will  be  determined  by 
the  area  ratio  of  the  inlet  duct  and  the  passage.  The  design  pressure 
ratio  for  the  blade  is  2.239  (21),  which  is  lower  than  the  pressure 
ratios  used  in  this  study. 
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TABLE  4 


COMPARISONS  TO  DESIRED  INLET  CONDITIONS 


*  ATEC  value 
••  ATEC  value 

A  trial  case  of  ANTEC  was  run  using  =.4  to  calculate  the 
necessary  input  values  of  T_  and  P„.  For  this  case,  the  flow  was  choked 
through  the  passage  and  the  inlet  values  do  not  drop  to  the  desired  FlOO 
values.  If  this  case  had  been  investigated  further,  choked  flow  should 
not  present  a  problem  for  this  particular  study.  If  matching  inflow 
conditions  were  obtained  from  ATEC,  a  valid  comparison  would  still  be 
possible.  Driver  and  Beran  have  found  in  their  research,  using  the 
rounded  trailing  edge  blade,  that  has  an  asymptotic  upper  limit  of 
.367  at  a  pressure  ratio  of  3.239.  (7:13) 

The  percentage  difference  is  presented  in  Table  4  as  a  point  of 
interest.  What  the  author  finds  curious  was  how  the  numerically 
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calculated  inlet  pressures  and  temperatures  equal  the  analytical  values. 
The  desired  inlet  conditions  were  based  on  the  values  of  M2  in  Table  2, 
and  the  values  of  M2  generated  by  the  codes  are  significantly  different. 
Since  emphasis  was  placed  on  matching  inflow  conditions,  the  difference 
between  the  values  of  Mj  was  not  considered  a  problem.  If  further  work 
is  done  on  this  subject,  the  FlOO  pressure  ratio  and  inlet  pressure 
could  be  replaced  by  blade  design  data.  (21)  Alternately,  the  blade 
could  be  replaced  with  the  FlOO  blade  geometry.  Eliminating  the  mixing 
of  geometries  and  conflicting  parameters  should  isolate  or  eliminate 
this  curiosity. 

III. 2  CFL  Criteria  and  Convergence  Histories 

The  Courant-Friedrichs-Lewy  (CFL)  condition  is  a  stability 
requirement  where  v,  called  the  Courant  number,  is  defined  by 

V  =  C  A  ^  /  Ax  (52) 

(1:75) 

All  ATEC  cases  used  a  CFL  value  of  .95.  However  to  successfully 
run  ANTEC  for  Cases  1  and  3,  the  time  step  had  to  be  adjusted  after  a 
number  of  iterations.  Following  start  up,  ANTEC  Case  1  yielded  a 
satisfactory  transient  solution  until  at  a  certain  time  (t  =  .047 
seconds)  the  numerical  solution  became  undefined,  thus  making  the  time 
step  suspect.  (5)  Therefore,  the  cases  were  then  run  in  intervals  of 

iterations,  as  demonstrated  in  Table  5.  If  the  solution  became 

undefined  in  an  interval,  the  CFL  number  was  decreased  and  the  code  was 

restarted  (see  Appendix  C  for  a  discussion  on  the  RESTART  file) .  After 
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adjusting  the  CFL  number,  the  code  again  would  progress .  Once  the 
solution  was  approaching  steady-state,  the  CFL  criteria  could  again  be 
raised.  A  similar  problem  was  encountered  by  Driver  and  Beran  in  their 
research  in  the  test  cases  where  the  contact  surface  was  in  the  vicinity 
of  the  trailing  edge.  (7:15)  For  this  reason,  it  was  originally  assumed 
that  the  numerical  solution  for  ANTEC  became  undefined  at  the  trailing 
edge.  However,  investigation  of  the  flow  showed  that  the  problem  was  at 
the  inlet.  Figure  3  shows  the  region  at  the  inlet  where  the  solution 
becomes  undefined.  Since  this  problem  was  only  seen  in  ANTEC  and  not 
ATEC,  the  inlet  boundary  condition  for  vibrational  energy  was  suspected 
as  being  the  cause. 


\  ‘ 

Figure  3:  Region  of  Undefined  Solution 


Table  5  summarizes  the  progression  of  CFL  numbers  for  the 
different  cases.  An  additional  high-temperature  case  is  listed  as  a 
point  of  interest.  Case  3  was  never  run  at  an  initial  CFL  value  of  .95 
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due  to  an  intermediate  case  with  Tg  =  2000°K  having  problems  after  1000 
iterations  at  .95.  Therefore,  it  is  not  known  if  Case  3  could  be 
started  at  a  higher  CFL  number  and  then  be  adjusted  as  with  Case  1 . 


TABLE  5 
CFL  CHANGES 


Case 

Number  of 

iterations 

CFL 

1 

1-3500 

.  95 

3501-4500 

.  90 

4501-8000 

.  95 

2 

1-8000 

.95 

3 

1-9000 

.  90 

9001-10000 

.95 

T2=3000 

1-1000 

.60 

1001-4000 

.50 

4001-7000 

.60 

7001-9000 

.70 

To  determine  if  steady-state  had  been  reached,  the  convergence 
history  for  each  case  was  recorded.  These  histories  are  shown  in 
Figures  4-9;  with  the  convergence  being  determined  by  the  normi 
calculated  as 


N(U)  =  100  {max  [abs(1  -  uI|/Uk*’)]}  (53) 

which  gave  a  percentage  difference  between  the  present  and  previous 
solutions  at  each  point  with  the  final  norm  value  being  the  maximum 
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difference  in  the  grid.  Therefore,  if  just  one  point  continually 
fluctuated  due  to  a  numerical  or  physical  phenomena,  the  convergence 
history  may  appear  not  to  reach  a  steady  state  as  seen  in  ANTEC  Case  3, 
Figure  9.  Had  the  nonequilibrium  Case  2,  Figure  7,  been  run  for  as  many 
iterations,  it  may  have  exhibited  similar  behavior.  These  oscillations 
are  consistent  with  Driver  and  Beran's  results  for  the  rounded  trailing 
edge  blade  used  in  this  study.  (7:10) 


Figure  4:  Convergence  History  for  ATEC  (Case  1) 


Figure  5:  Convergence  History  for  ANTEC  (Case  1) 
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Figure  9;  Convergence  History  for  ANTEC  (Case  3) 
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Except  for  a  spike  in  the  ANTEC  results,  both  codes  had  identical 


convergence  histories  for  Case  1  (Figures  4  and  5) .  For  Cases  2  and  3, 
the  norms  from  ATEC  reached  a  'smoother'  steady-state  solution  than 
those  from  ANTEC  (Figures  6-9) .  In  Figure  9,  the  oscillations  appear  to 
begin  where  the  CFL  was  raised  to  .95.  However,  there  is  no  such 
correlation  between  CFTi  and  the  norms  in  Figure  5.  Therefore,  more  test 
cases  would  have  to  be  run  before  an  assertion  of  a  connection  between 
CFL  number  and  solution  oscillations  was  made. 


_ Blade  Temperature  Profiles 

The  trends  observed  for  the  comparisons  of  the  blade  temperature 
profiles  did  not  correspond  to  those  of  a  blunt-nosed  body  in  hypersonic 
flight  as  previously  mentioned  in  Section  I.l.  However,  for  the  flow 
conditions  present,  the  results  were  realistic. 

Figures  10-12  are  the  blade  temperature  profiles  for  the  three 
cases  investigated.  The  stagnation  temperature  on  the  lower  blade 
surface,  known  as  the  pressure  surface,  was  as  expected  and  corresponded 
closely  to  the  inlet  temperature.  For  Case  2,  the  results  from  ANTEC 
show  the  stagnation  temperature  was  1673.  l^^K,  or  0.54%  different  than 
the  inlet  temperature;  and  for  Case  3,  it  was  2484. 37°k,  0.63%  different 
(Figures  11  and  12) . 

At  the  blade  tip,  there  was  a  very  close  correspondence  between 
the  blade  temperatures  predicted  by  the  two  codes,  which  was  also 
expected.  This  part  of  the  blade  was  in  equilibrium  as  shown  by  the 
overlapping  of  the  blade  temperature  and  from  ANTEC.  Here,  the 
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internal  energies  of  the  flow  produced  by  the  two  codes  were  equal. 


With  Bin,  defined  as 


e.n.  =  cj 


the  following  relationship  was  observed 


(54) 


ne 

=  (CyT)a 

(55) 

For  Cases 

2  and  3,  as 

the  flow 

moved 

from 

this  reference 

point. 

expanding 

and  cooling 

over  the 

upper 

(or 

suction)  surface 

,  the 

temperature  decrease  behaved  differently  for  the 

two  codes  based 

on  the 

values  of 

Yeff  Ya-  Yeff/a 

(Equation 

14)  can 

also 

be  defined  as 

R 

7eff/a  ^  f.  (55) 

i-v 

With  Ya  ^  constant,  Cy  must  also  be  constant  (a  calorically  perfect 
gas)  .  In  ANTEC,  neither  Yeff  ^  assumed  constant,  but  were 

temperature  dependent  as  demonstrated  in  Equations  14-17.  The  value  of 
Cy  denotes  the  capacity  of  the  gas  to  hold  energy.  As  the  temperature 
decreased,  ATEC  retained  the  same  capacity  to  store  energy.  However,  in 
the  nonequilibrium  version,  Yeff  increasing  due  to  Cy  decreasing; 

capacity  to  store  energy  was  lost,  resulting  in  higher  temperatures. 

Another  way  to  look  at  this  is  to  assume  the  flows  in  both  code 

versions  were  transferring  internal  energy  to  kinetic  energy  at  the  same 

rate  along  the  blade,  retaining  internal  energies  that  were 

approximately  equal  to  each  other.  At  any  arbitrary  point  along  the 

blade  (CyT)^  =  (CyT)g.  Again  observing  the  suction  surface,  the  value  of 

C«  was  constant,  but  C,  had  decreased  due  to  the  cooling  in  the 
’a  ne 
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expansion.  Clearly  for  the  energies  to  be  equal,  T„g  must  be  greater 
than  Tg. 
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Figure  10:  Blade  Temperature  Profile  Case  1 
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Figure  11:  Blade  Temperature  Profile  Case  2 
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Figure  12:  Blade  Temperature  Profile  Case  3 


The  blade  temperatures  for  the  two  methods  at  the  low-temperature 
case  were  as  expected,  very  close  to  equal.  However,  the  low- 
temperature  case  posed  a  problem  to  the  nonequilibrium  model.  The 
characteristic  time  of  the  problem,  was  smaller  than  the  relaxation 
time,  T.  Physically,  for  a  cold,  low-pressure  flow,  relaxation  will  be 
relatively  long.  For  this  case,  at  8000  iterations,  was  on  the 
order  of  lO"^  and  t  was  on  the  order  of  10^.  Therefore,  the  flow  had 
not  reached  equilibrium.  Clearly,  the  solution  would  be  to  employ  local 
time  stepping.  By  Equation  8,  as  temperature  increases,  t  will 
decrease.  Since  t  for  Case  2  was  on  the  order  of  10~'°,  no  time  scale 
problem  existed  for  the  high-temperature  cases  and  the  nonequilibrium 
model  was  able  to  reach  a  physical  solution. 
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:iy  £.utje  lemueraturtia 

After  Studying  the  overall  blade  temperature  profile,  the  trailing 
edge  was  of  particular  interest  due  to  the  complex  flow  structure  in 
this  area. 

On  the  pressure  surface,  the  temperature  was  relatively  constant 
until  the  trailing  edge  (Figures  10-12) .  These  blade  profiles  showed  a 
rapid  expansion  and  decrease  in  temperature  followed  by  an  immediate 
compression  or  shoc)c.  The  suction  surface  also  had  a  rapid  compression 
at  the  trailing  edge.  The  changes  in  temperature  before  this 
compression  were  due  to  the  blade  and  cascade  geometry.  The  shock 
created  at  the  trailing  edge  of  the  pressure  surface  was  impacting  the 
suction  surface,  causing  a  slight  compression.  An  expansion  was  caused 
by  the  rounded  trailing  edge.  The  rapid  compressions  experienced  by 
both  surfaces  were  in  order  to  meet  the  Kutta  condition,  which  states 
that  pressures  from  the  suction  and  pressure  surfaces  must  be  equal  at 
the  trailing  edge  stagnation  point. 

Figures  13  and  14  highlight  the  trailing  edge  geometry  and  show 
the  differences  between  the  two  code  versions  in  the  vicinity  of  the 
trailing  edge  for  Case  3.  For  this  case,  the  maximum  temperature 
differences  over  the  entire  grid  were  -741'’K  and  539°K.  These  were  at 
the  trailing  edge  as  demonstrated  by  the  extreme  light  and  dark  shaded 
regions  in  Figure  13  and  were  not  evident  on  the  blade.  Although 
weakened.  Figure  14  also  demonstrates  the  pressure  surface  trailing  edge 
shock  impacting  the  suction  surface.  The  shape  of  the  trailing  edge 
appears  different  in  these  figures  due  to  the  different  scalings. 
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III.5 _ Gamma  Comparisons 

Figures  16  and  17  show  how  gradually  changes  along  the  blade 
(at  the  ghost  points)  with  temperature,  and  it  can  be  seen  that  the 
value  of  everywhere  on  the  blade  was  higher  than  y^.  Therefore,  even 
through  the  shoc)c,  where  temperatures  sharply  increase,  the  increase  was 
not  enough  to  significantly  increase  the  capacity,  Cy,  of  the 
nonequilibrium  flow  and  cause  lower  temperatures  than  those  predicted  by 
ATEC. 

Again,  the  low-temperature  case  (Case  1)  showed  difficulty  due  to 
the  flow  not  being  able  to  reach  equilibrium.  An  accurate  value  of 
could  not  be  calculated.  The  expected  trend  as  the  flow  expanded  would 
be  for  Ygff  to  approach  1.4.  Figure  15  shOKS  this  is  not  the  case. 

The  blade  profiles  of  Ygfj  for  all  cases  closely  resembles  the  Ty,;, 
curves  in  Figures  10-12,  but  inverted.  For  ease  of  comparison.  Figure 
18  is  the  vibrational  temperature  for  Case  3.  In  ANTEC,  is 

extracted  from  e^it,  (Equation  7)  and  6^^,  is  dependant  on  (Equation 

53)  .  Knowing  the  dependence  of  Yg{f  on  a  relationship  between  Yett 

and  Tyj(,  can  be  deduced  showing  that  as  Yett  increases,  Ty^,  decreases. 
This  trend  is  indeed  shown  between  Tyig  and  Ye«  the  blade.  Figures  17 
and  18.  The  decoupling  of  vibrational  energy  does  indeed  'correct'  the 
value  of  Yett  i  • to  an  appropriate  value  for  the  flow. 
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Figure  15:  Gamma  Comparison  for  Case  1 


Figure  16:  Gamma  Comparison  for  Case  2 


X  Along  the  Blade 


III. 6  Computational  Time  Comparisons 


As  important  as  accurate  models  are,  there  is  usual  Iv  =  .jorialty 
paid  with  increased  CPU  time.  Most  of  the  runs  for  the  cases 
investigated  were  run  in  the  "background"  system  mode.  Although  an 
efficient  way  to  run  a  large  job,  1000  iterations  could  take  anywhere 
from  1  to  2  1/2  hours  to  complete,  depending  on  system  usage. 

ATEC,  and  hence  its  modification,  ANTEC,  were  optimized  to  run  on 
a  CRAY-XMP;  however,  all  computations  for  this  effort  were  performed  on 
Stardent  ST2000.  Therefore  parallelization  on  the  Stardent  was  not  an 
efficient  use  of  the  processors  and  the  runs  were  completed  using 
vectorization  only. 


TABLE  6 

COMPUTER  TIME  COMPARISONS 
CPU  SECCNDS/ITERATION 


Version 

No 

Optimize 

Vectorize 

Vectorize  and 
Parallelize 

ATEC 

7 . 53 

0.84 

1 . 99 

ANTEC 

11.79 

1.36 

5.93 

Interactive  runs  were  made  to  compare  the  run  times  of  ATEC  with 
those  of  ANTEC  at  different  levels  of  optimizations.  Table  6  summarizes 
the  average  CPU  time  per  iteration  for  the  two  versions.  Both  were 
compiled  and  run  using:  1.  no  optimization;  2.  vectorization  only;  ana 
3.  parallelization  and  vectorization.  With  no  optimization,  the  codes 
were  only  run  for  100  iterations.  The  codes  comipiled  this  way  were 
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never  used  to  arrive  at  a  steady-state  solution  of  the  gasdynarrdc 
equations.  The  other  test  runs  were  performed  with  1000  iterations  and 
are  more  representative  of  actual  'runs.’  The  time  required  to  output 
the  data  files  after  completion  of  the  iterations  can  considered 
constant;  therefore,  these  average  times  per  iteration  should  decrease 
slightly  as  the  number  of  iterations  is  increased. 
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This  comparative  study  between  ATEC,  a  two-dimensional,  explicit 
Euler  code,  and  a  version  modified  to  account  for  vibrational  energy, 
ANTEC,  has  shown  the  impact  of  y  in  computing  transonic  flow  through  a 
turbine  cascade.  The  assumption  of  a  calorically  perfect  gas  may 
simplify  the  computations;  however,  even  when  corrected  for  high 
temperatures,  a  constant  value  of  y  for  the  entire  flowfield  is 
inaccurate.  It  was  the  purpose  of  this  investigation  to  determine  the 
limitations  of  the  assumptions  and  the  models.  For  Case  3,  the  extreme 
temperature  differences  were  -741®K  and  539®K  and  were  located  at  the 
trailing  edge.  These  values  may  not  appear  as  significant  as  the 
5,100°K  temperature  difference  for  the  hypersonic  blunt-nose  body 
mentioned  in  Section  I.l;  however,  such  extreme  differences  were  not 
anticipated. 

The  use  of  ANTEC  for  the  high-temperature  cases  resulted  in  higher 
blade  temperatures  over  the  suction  surface.  With  steady,  equilibrium 
inflow  conditions,  the  tip  of  the  blade  was  also  in  equilibrium  as 
demonstrated  by  T  =  T^ib  in  this  area.  Areas  of  nonequilibrium  were  seen 
toward  the  trailing  edge  on  the  suction  surface  and  at  the  trailing 
edge . 

Even  at  low  temperatures,  where  the  nonequilibrium  model  showed 
difficulty  due  to  a  large  relaxation  time,  there  were  small  differences. 
However,  these  differences  could  be  considered  minor  depending  on  the 


accuracy  of  the  estimation  desired.  Having  based  the  algorithm  of  the 
Euler  equations,  other  important  aspects  of  transonic  flow  through  a 
turbine  cascade  such  as  heat  transfer  are  ignored,  making  the  code  a 
tool  for  estimation.  Incorporation  of  the  thermal  nonequilibrium  model 
allows  the  designer  to  simulate  higher  temperature  flows  with  a  greatly 
improved  level  of  accuracy. 

IV. 2  Recommendations  for  Further  Study 

1.  Eliminate  the  mixture  of  blade  and  design  parameters  and  vary 
the  inlet  temperatures.  This  would  possibly  identify  any  peculiarities 
in  this  investigation  caused  by  using  FlOO  inlet  parameters  with  the 
NASA  test  blade  geometry. 

2.  Modify  ANTEC  to  use  a  steady-state  ATEC  RESTART  file  as  the 
initial  conditions  to  help  alleviate  the  computer  time  required  for  the 
nonequilibrium  model. 

3.  Incorporate  the  nonequilibrium  model  without  the  assumption  of 
air  being  a  single  species. 

4.  Conduct  a  parametric  study  varying  input  parameters  discussed 
in  Appendix  C  that  were  held  constant  for  this  investigation. 

5.  To  simulate  equilibrium  flows,  modify  ATEC  to  have  a  'look-up' 
table  for  Yg  and  compare  the  results  and  computer  times  with  the 
nonequilibrium  version. 
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To  gain  an  understanding  of  second-order  accurate  shock  capturing 
schemes  and  the  nonequilibrium  model,  unsteady  flow  in  a  shock  tube  was 
investigated.  The  shock  tube  is  a  valuable  learning  tool  since  an 
analytic  solution  can  be  found  and  much  is  available  in  the  literature 
exists  for  comparisons. 

Gas  in  the  left  side  of  the  tube  was  initially  set  to  high 
pressure  and  density  with  the  other  side  being  of  lower  pressure  and 
density.  Ideally,  a  shock  should  be  captured  by  no  more  than  one  node, 
since  this  would  more  closely  model  the  instantaneous  changes  caused  by 
a  shock.  Therefore,  capturing  the  shock  with  a  minimal  number  of  nodes 
is  desirable.  Also,  a  shock  tube  not  only  has  the  shock  discontinuity, 
but  changes  in  flow  parameters  due  to  the  contact  surface  and  expansion 
fan  as  well. 

A  comparison  of  the  nonequilibrium  code  versus  an  ideal  gas  code 
was  made  (Figure  19)  .  The  nonequilibrium  results  nearly  match  those 
reported  by  Moran.  (16,17) 
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FIGURE  19:  Shock  Tube  Results 
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The  TVD  scheme  used  in  Phases  One  and  Two  was  different  from  the 
one  applied  in  ATEC.  The  shock  tube  used  Marten's  ULTIC.  (12)  The 
differences  between  the  two  schemes  are  now  described. 

There  is  no  change  to  the  form  of  the  finite  difference  scheme  or 
numerical  flux  (equation  numbers  are  repeated  for  identical  formulas) 


U""'  =  U" 


7  n 

I  (f 


H-l/2 


-  t  Im)  + 


AtS 


(32) 


f 

\*\a 


F(Up  +  F(Uj,,)  + 


>1/2 


\*\a 


k=l 


(33) 


is  still  the  eigenvector  matrix  (Equation 

numerical  viscosity,  B  t  is  defined  as 

>1/2 


22)  but  the  effective 
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(57) 


\*\a 


=  (9, 


Q(v^) 


a, 


■j+1/2 


with  the  same  values  of  v  ,  a^|,2'  Q 


''2  “  ^  ^1*1/2  \\a 


(35) 


''1  “  ^  ^1+1/2 


(36) 


“j»i/2  “  [^i+w] 


(37) 


Ah/2U  =  U,,,  -  Uj 


(38) 


Q(Z)  =  "^  +  e  |z|  <  2e 


Q(Z)  =  lz|  for  |z|  >  2e 


(40) 


but  a  becomes 


^1/2(2)  =1  (1  -  Q(z) ) 


(58) 


e  is  still  an  input  parameter.  In  the  shock  tube  code  a  value  of  1/2 
was  used  due  to  lesser  values  causing  an  undefined  solution  (a  value  of 
0  was  used  in  both  ATEC  and  ANTEC) .  This  sensitivity  was  also  seen  by 
Moran.  (16) 

The  linearly  degenerate  fields,  a*' =  u,  are  described  first,  gl*  is 
determined  from 
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(59) 


-K  — K  ^  — K 

g,  =  Qj  +  9j 


9]  ®f+l/2  (IQj+mI  /  9j.i;2  ^j*1/2  (60) 

9]  ~  ^f+1/2  rr'3x[  0,min  (Sj^],2  Oj.1,7  ®j-i,7' ®i«i/2  l^j-i/2l*  1  (61) 

gil/2  =  |[Q(v^)  -  (v'!yK,a  (62) 

with 

S|ii/2  =  signd,  ai,/^)  (42) 

~  sign  ( 1 , Qj^,^)  (63) 


where  'sign'  operator  assigns  the  sign  (+  or  -)  of  the  second  value  in 
parenthesis  to  the  first  value. 


Now  the  treatment  of  the  nonlinear  fields,  a*' =  u  ±  C,  is  described.  gj'  is 
determined  from 


0  is  defined  by 


l%17  ■ 


l®j  I/2I 


•j-1/21 


(64) 


-k 

9i  = 


-k 


(65) 


is  unchanged  from  the  Harten-Yee  formulation  used  in  ATEC  and  ANTEC 
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^9j+1/2  9j-1.'2^ 

\V2  =  J “'i/2  ^  0  (^^a) 

“)+l/2 

I 

ylv2  =  0  %xa  =  0  (4  4b) 

The  values  to  be  used  for  the  eigenvalues  and  eigenvector  matrices 
for  the  equilibrium  code  (Phase  One)  can  be  found  in  references  4  and 
12,  and  in  references  15,  16  and  17  for  the  nonequilibriumi  code  (Phase 

Two)  . 
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Due  to  the  modularity  of  ATEC,  incorporation  ol  the  nonequilibrium 
model  into  ATEC  was  very  straightforward.  The  largest  task  was  changing 
the  Jacobian  matrices  and  their  inverses  to  account  for  the  added 
vibrational  energy  equation;  systems  of  4x4  matrices  were  extended  to 
systems  to  5x5  matrices. 

Along  with  p,  pu,  pv  and  E,,  the  vibrational  energy  had  to  be 

calculated  at  cell  centers,  and  Roe  averaged  at  cell  interfaces.  From* 
®vib'  "'"vib  could  be  calculated  using  Equation  7.  was  not  necessary  to 

obtain  a  steady-state  solution,  but  it  was  necessary  to  determine 
whether  or  not  the  solution  was  in  equilibrium  once  steady  state  was 
attained. 

Calculation  and  application  of  the  source  term  was  the  significant 
difference  between  ANTEC  and  ATEC.  Therefore,  Equations  2,  6,  8,  9,  10, 
and  13  were  incorporated  into  ANTEC  along  with  the  modified  equations 
for  pressure,  total  energy  and  speed  of  sound  (Equations  12,  3,  and  28) . 

To  verify  proper  implementation  of  the  nonequilibrium  model,  the 
modified  version  was  run  with  the  source  term  set  to  zero  (conservation 


5rational  energy) 


identical  to  those 


unmodified  version  as  seen  in  Figure  20. 


ATEC  and  ANTEC  are  general  codes  written  to  run  for  a  variety  of 
flow  conditions  and  grid  sizes,  as  dilineated  in  the  CASIN  file. 

This  file  first  defines  and  T_.  These  are  derived  using 

Equations  47  and  48.  To  obtain  an  accurate  comparison  between  ATEC  and 
ANTEC,  it  was  important  to  match  the  inlet  conditions  for  each  test 
case . 

Next  CASIN,  lists  the  flow  angle.  This  parameter  was  not  altered 
for  any  of  the  cases. 

The  plenum  pressure,  obtained  from  calculations  based  on  the 
desired  inlet  pressure  and  stage  pressure  ratio,  is  then  input. 

The  following  two  fields  are  used  to  start  the  codes  running  from 
zero  or  to  use  a  RESTART  file  and  determine  how  many  iterations  to  run. 
If  a  1  is  entered,  the  computational  domain  will  be  initialized  and 
computations  will  begin.  If  any  number  other  than  1  is  listed,  a 
RESTART  file  will  be  used.  The  number  of  iterations  specified  will  be 
run  and  the  U  vector  for  all  grid  points  will  be  written  to  RESTART. 
With  the  restart  capability,  the  code  can  be  stopped,  the  development  of 
the  flow  field  can  be  examined,  and  then  the  code  can  be  started  again. 
This  feature  is  also  convenient  when  a  steady-state  solution  is  desired, 
yet  how  many  iterations  are  necessary  to  reach  this  condition  is  not 
known.  If,  for  example,  5000  iterations  are  accomplished  and  the  code 
i?  stopped  and  the  solution  is  not  yet  at  steady- state;  then  the  CASIN 
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is  easily  edited  for  more  iterations  using  the  solution  at  5000. 
Without  the  restart  capability,  much  computer  time  would  be  wasted. 

Next,  the  CFL  condition  is  input.  For  ATEC  this  was  always  .95. 
However,  this  parameter  had  to  be  adjusted  for  ANTEC  which  was  discussed 
(  in  Section  III. 2. 

The  code  can  run  either  a  first  order  or  second  order  solution. 
With  and  input  value  of  1,  the  code  will  run  a  first  order  Roe  scheme; 
with  a  value  of  2,  the  code  will  run  the  second  order  TVD  scheme  based 
on  ULTIC.  The  value  2  was  always  used  in  this  investigation. 

The  next  seven  fields  deal  with  the  grid  by  defining  the  maximum 
number  of  nodes  in  each  direction,  the  location  of  the  diaphragm,  and 
the  positions  of  the  blade  trailing  edge  and  the  periodic  boundaries. 
These  parameters  were  never  varied. 

As  mentioned  in  the  section  dealing  with  the  development  of  the 
TVD  scheme,  e  and  (o  are  input  parameters.  These  also  were  never 
varied.  The  value  of  e  was  zero  for  both  the  linearly  degenerate  and 

genuinely  nonlinear  fields  and  (0  was  equal  to  2,  as  used  by  Yee,  for 

the  linearly  degenerate  fields  and  zero  for  the  nonlinear  ones.  (23:28) 
Lastly  the  grid  file  is  defined.  Again,  the  same  grid  was  used 
for  all  cases  in  this  study. 

As  a  final  note,  if  the  original  version  of  ATEC  is  to  be  run  at 
high  temperatures,  using  a  temperature  corrected  value  of  rather  than 
the  ideal  gas  value  of  1.4,  y^  should  be  added  as  an  input  rather  than 

'hard  wired'  into  the  code.  This  would  save  having  to  edit  and 

recompile  for  different  temperatures. 

The  following  are  sample  CASIN  files  for  the  two  versions 
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‘Jample  CASIN  for  ANTEC 


320.24e4  2798.41 
49.2 
58.96e4 
2  10003 
0.95 
2 

177  20  0 
79  99  40  138 
0.0  0.0  0.0  0.0 
2.0  2.0  2.0  0.0 
blntgrd 

Sample  CASIN  for 
311.30e4  2714.37 
49.2 

58.965e4 
2  7001 
0.95 
2 

177  20  0 
79  99  40  138 
0.0  0.0  0.0  0.0 
0.0  2.0  0.0  2.0 
blntgrd 


(PtINF,  TtINF) 

(BETA2) 

(P3) 

(CEL,  LOCAL) 

(METHOD;  l=Roe,  2=ULT1C) 
(IMAX,  JMAX,  IDIS) 

(IPERl,  IPER2,  ITEl,  ITE2) 
1.0  (EPSILON  1,  2,  3,  4  and  5) 

1.0  (OMEGA  1,  2,  3,  4  and  5) 

(GRID  FILE) 

the  original  ATEC 

(PtINF,  TtINF) 

(BETA2) 

(P3) 

(CFL,  LOCAL) 

(METHOD:  l=Roe,  2=ULT1C) 
(IMAX,  JMAX,  IDIS) 

(IPERl,  IPER2,  ITEl,  ITE2) 
(EPSILON  1,  2,  3  and  4) 
(OMEGA  1,  2,  3  and  4) 

(GRID  FILE) 
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